Importance

Bioturbation is the disturbance of sediments by biogenic reworking, and biodeposition is the deposition of biogenic material in the sediment.

These two processes are important in the soil and benthos dynamics because they increase ground surface area and particle mixing, which in turns enable and enhance ground oxygenation, gas interchange and nutrient recycling.

Crab bioturbation has been found to influence the vegetation structure, plant productivity and the bacterial community composition in the sediment.

Background

Traditionally, bioturbation in burrowing crabs is studied by analyzing burrow casts or the amount of loose sediment in the surface around a burrow. Both methods assume that bioturbation and biodeposition are a function of sediment removed and deposited elsewhere, and sediment compactation.

The burrow casts method uses the total volume of the burrow chamber and its surface area as a surrogate of bioturbation. Estimates of bioturbation rates (i.e. volume per unit time) using this approach are difficult and prone to speculation as the absolute volume of the chamber represents the sediment moved and compacted over an undertermined number of days. Given its invasive and costly nature, this method is only applied to few burrows. Thus, bioturbation rates are calculated by extrapolating estimates based on the total number of burrows and the estimate volume. Variable burrow diameter and architecture can introduce additional uncertainties in this method. The study of burrow casts is worth persuing, and it can inform and complement other methods, but by itself is likely to produce bioturbation and biodeposition estimates with large uncertainty intervals.

The amount of loose sediment is a direct proxy of the amount of sediment reworked by crabs. By using this method it is possible to obtain precise values of the amount of mass moved by crabs. This method requires frequent sampling, so loose sediment created by crab activities, can be identified and collected before is eroded or compacted. Collecting loose sediment in the surface around a burrow openning requires little equipment and intervation, so it can be done over relatively large areas. Collecting the loose sediment deposited inside the burrows is more complicated (as far as I know it has not been done). Loose sediment sampling from inside the burrow is likely to affect the structurtal integrity of the chamber, and produce additional loose sediment which can increase sampling error. For this reason, some researchers (e.g. Wang et al 2010) have used artificial burrows (made of PVC) to measure the amount of loose sediment deposited inside burrows by crabs.

Note: Using burrows made of PVC to quantify loose sediment deposited inside burrows has its own issues which I am ignoring at this stage.

Collecting loose sediment is a better solution for assesing crab bioturbation. Sediment collected can also be used to answer additional questions about the consequences of bioturbation in the sediment properties. While this method require little equipment during sampling, it can be very time consuming when applied over large areas or areas with high burrow density.

The general assumption when studying crab bioturbation is that making a burrow requires excavating, compacting and/or depositing sediment in the burrow or elsewhere. Thus, the number of burrows and their size are believed to be related to bioturbation and biodeposition rates. In addition, crab feeding activity will also contribute to bioturbation and biodeposition.

In this paper, the method described in Herrera et al (2020)) is used to estimate bioturbation and evaluate potential process underlying the observed bioturbation patterns.

Aims

  1. Assess the bioturbation and biodeposition activity of burrowing crabs.

    • Q: How much sediment matter (estimated as volume) do crabs move per unit time?

  1. Study the potential process explaining bioturbation patterns.

    • Q: Do the number of burrows explain the observed change in volume?

    • Q: Do the distribution of burrows in the quadrat is associated with areas of major change in volume?

    • Q: Do the burrows distribution correlates with the local geometry of the sediment surface? This is, burrows position explain the observed changes in verticality, surface variation and other geometric features.

Methods

We proposed (Herrera et al 2020) a new method using computer vision and photogrammetry to assess the amount of volume change due to burrowing and feeding activities by crabs.

Our method is a fundamental improvement to the estimation of crab bioturbation as it relates sediment turnover (i.e. moved matter) to observed volume changes in the sediment surface per unit time.

This study was performed in a mudflat next to the Ross River in Townsville, Queensland, Australia. The area was selected considering easy access and the proven presence of fiddler crab from species Tubuca polita and Tubuca coarctata. The tidal regime in this area is well known and inundation can be predicted with precision.

Five permanent quadrats were installed in sediment areas exclusively or predominantly dominated by Tubuca polita individuals. Tubuca polita individuals are primarily active in the sediment surface during/following spring tides. To avoid confounding bioturbation with other natural process that can modify or moved matter, such as tidal erosion, the permanent quadrats were installed out of the area of tidal influence, this means, quadrats were not flooded.

All quadrats covered a 1 $m^{2}$ (1 x 1 m). Vertices from each quadrat were marked with Ground Control Points (GCPs). During the sampling time only two GCPs were lost. Replacements were installed inmediately after a lost was identified using a plastic reference quadrat.

Quadrats were followed during 14 days, and sampling occurred simultaneosly for all quadrats, same day with minutes of difference. First sampling event occurs at day zero (DO), with succesive samplings after one day (D1), three days (D3), seven days (D7), ten days (D10), and fourteen days (D14). Each sample event consisted in recording a video perpendicular to the sediment following the method described in Herrera et al 2020.

Unfortunately, due to unforeseen equipment failure samples for D3, D7 and D10 were not analyzed due to file corruption or due to modified recording settings which made results incompatible. Thus, this study only presents results for samples D0, D1 and D14.

Important: The issues associated to days D3, D7 and D10 are largely my responsability because I did not check the videos just after recording, instead I did just before sampling D14 which allowed me to identify the problem and fixed it before D14.

For each quadrat and sample event three products were created using Structure from Motion following Herrera et al (2020) methodology: point cloud data, textured and coloured mesh model and coloured orthomosaics.

The difference in volume for each quadrat and sample event was calculated by comparing point cloud data between pair of sample events.

Point Cloud Data

Point Cloud Data refers to coordinate data (x, y, z) created from Photogrammetry and other techniques. Point Cloud Data is the essential data type from which many 3D modelling and analysis tools based upon. This type of data represents a point in space based on its coordinates x, y, and z in a 3D space.

As any other data set in higher dimensions (> 2D) statistical techniques can be used to reduce dimensionality. PCA is used in during the analysis of point cloud data in many computer vision applications (e.g. find features such as the edges of a geometry, or define the slope in an particular area).

Below there is a small representation of a point cloud data set. We could calculate a PCA for the entire dataset, and it will result in a new representation of the data in a 2D space. It is very likely that in such PCA the two principal components will be associated to axes "x" and "y", as much of the variation occurs in these two axes (i.e. there is not much changes in "z" or elevation). Thus, eigenvalues for "x" and "y" will be larger than "y" eigenvalues.

Tip: You can use your mouse to rotate the view and zoom in and out.

Alternatively we could recursively calculate PCAs in small windows including just a small subset of the data. Thus, we will obtain many PCAs for small sections of the data. Areas where most variation occurs in the plane "xy" will have show larger eigenvalues in "x" and "y". However, if the window is small enough, we could detect areas where change in elevation ("z") are more important than changes in the plane "xy". Therefore, "z" eigen value will be larger.

Based on this principle, we can reduce our point data dimensionality to 2D (i.e. an image or raster). Based on linear combinations of the eigenvalues, computer scientists have came with multiple metrics that describe geometric features of point cloud data. Some of these are: PCA1, PCA2, Planarity, Spherecity, Surface Variation, Verticality and 3rd Eigenvalue.

These features can be relevant when describing biological processes, such as crab activities.

How do these features look like?

The geometric features described above are presented in the compound figure below for one video sample.

Note: A video sample represents one video scan for one quadrat at a given time.

The first plot (top-left) shows the elevation (or height in cm) for this quadrat. Elevation represents the value of the "z" coordinate. The second plot (top-middle) shows roughness, which is a measure of the distance (cm) between a specific point and the plane made by the neighbors points inside a specific sphere window (radius = 1 cm). Essentially high roughness values indicate uneven or irrfular surfaces, i.e. points which are distance from the plane created by its neighbors points. All other plots represent geometrics features that were calculated after calculating PCAs for a sphere window (radius = 1 cm). Features are arithmetic combinations from eigenvalues and eigenvectors (as per Hackel et al 2016).

Some geometric features are not very informative in this particular tasks. However, some can provide interesting insgights into the local shape/geometry in the surface. For instance, verticality (top-right) can be interpreted as changes in the sediment slope. Planarity, Surface Variation and Sphericity are also good descriptors of the sediment roughness.

#collapse-hide

options(repr.plot.width=15, repr.plot.height=15)

# Visualize burrows in quadrat and change in Digital Surface Models
## Load clean data
features_list <- readRDS(here::here('data/clean', 'features_list.rds'))

A_D0 <- features_list[[1]]

features <- names(A_D0)[c(4:length(names(A_D0)))]
plots <- list()
for (feat in features){
# Geometry feature plot
plots[[feat]] <- ggplot(A_D0, aes_string("x", "y", fill=feat)) +
  geom_tile() +
  scale_fill_viridis_c(name = paste0(feat)) +
  coord_equal(xlim = c(0, 100), ylim = c(0, 100)) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  labs(x = feat) +
  theme(# plot.margin = margin(5.5, 40, 5.5, 5.5), 
    legend.position = "right",
#     legend.key.size = unit(0.8, "cm"),
    legend.title = element_blank(),
    axis.title.y = element_blank(),
    # legend.justification = c("right", "top"),
    panel.background = element_blank(),
    axis.line = element_line(color = "black", size = .5),
    panel.grid.major = element_line(color = "gray", size = 0.25),
    legend.key = element_rect(fill = "white", color = NA))

}

plot_grid(plotlist = plots)

Textured and coloured mesh model

A second type of product that can be obtained from this workflow are textured and coloured mesh. In the current workflow, these are only used for educational purposes (e.g. visualization for presentation). The animation below shows one example of this product.

Note: This animation uses a simplified version of the data. Thus, a mesh with higher resolution could be created.

Coloured orthomosaics

Coloured orthomosaics, a 2D image or raster, are created by combaning multiple information (i.e. pixels) from the video scan. This images were used to count the number of burrows within each quadrat at each time.

Analysis and results

Aim 1: Assess the bioturbation and biodeposition activity of burrowing crabs

  • Question: How much sediment matter (estimated as volume) do crabs move per unit time?

To answer this question point cloud data was scaled and ortorectified. Then each possible pair of cloud data for a given quadrat (i.e. D0 vs D1, D0 vs D14, D1 vs D14) were aligned and superimposed using GCPs as coordinate references. The distance between the cloud pair was computed. Three comparisons per quadrat were possible: D0 was compared to D1 and D14, and D1 was compared to D14. These comparisons resulted in volume change estimates for three time periods: one, thirteen and fourteen days.

Important: the issues with cameras described before really messed my original plan of having more comparisons. 

#collapse-hide
options(repr.plot.width=10, repr.plot.height=10)
# Visualize burrows in quadrat and change in Digital Surface Models
## Load clean data
biotur <- read.csv(here::here('data/raw', 'bioturbation.csv'), header = TRUE, sep = ",",
                   colClasses = c("factor", "factor", "numeric", "integer", "numeric", "numeric", 
                                  "numeric", "numeric", "numeric", "factor", "integer", 
                                  "numeric", "numeric", "numeric"))

biotur$quadrat <- as.factor(as.character(stringr::str_sub(biotur$Cloud.1, start = 5, end = 5)))                
biotur$time1 <- as.factor(as.character(stringr::str_sub(biotur$Cloud.1, -7, -1)))
biotur$time2 <- as.factor(as.character(stringr::str_sub(biotur$Cloud.2, -7, -1)))
biotur$time1 <- sapply(biotur$time1, function(x) gsub("23NOV19", "D0", x), USE.NAMES = TRUE)
biotur$time1 <- sapply(biotur$time1, function(x) gsub("24NOV19", "D1", x), USE.NAMES = TRUE)
biotur$time1 <- sapply(biotur$time1, function(x) gsub("07DEC19", "D14", x), USE.NAMES = TRUE)
biotur$time2 <- sapply(biotur$time2, function(x) gsub("23NOV19", "D0", x), USE.NAMES = TRUE)
biotur$time2 <- sapply(biotur$time2, function(x) gsub("24NOV19", "D1", x), USE.NAMES = TRUE)
biotur$time2 <- sapply(biotur$time2, function(x) gsub("07DEC19", "D14", x), USE.NAMES = TRUE)
biotur$Bio_rate <- as.numeric(biotur$Bio_rate)

biotur$Time <- as.factor(as.integer(biotur$Time))

library(reshape2)
keep_names <- names(biotur)
keep_names <- keep_names[c(-5, -6)]
# keep_names

biotur1 <- melt(biotur, id.vars = keep_names)
biotur1$value <- ifelse(biotur1$variable == "Removed.volume", biotur1$value * -1, biotur1$value)


p1 <- ggplot(biotur1, aes(x = Time, y = value, fill=variable)) +
  geom_bar(stat = "identity") +
  facet_grid(~quadrat) +
  scale_fill_manual("Volume", values = c("#848484", "#000000"), labels = c("Added", "Removed")) +
  scale_y_continuous(breaks = seq(-3000, 9000, 1500)) +
  labs(x = 'Time lapsed (days)', y = 'Change of volume in the sediment (cm^3)') +
  theme(
    axis.text = element_text(size = 15),
    axis.title = element_text(size = 18),
    # legend.position = "none",
    panel.border = element_blank(),
    panel.background = element_blank(),
    panel.grid.minor = element_blank(),
    panel.grid.major = element_blank(),
    axis.line = element_line(color = "black", size = .5),
    # legend.key = element_rect(fill = "white", color = NA),
    axis.text.x = element_text(hjust = 1),
    strip.placement = "inside",
    strip.text = element_text(size = 18),
    strip.background = element_blank(),
    legend.text = element_text(size = 18),
    legend.title = element_text(size = 20))

p1


#collapse-hide
options(repr.plot.width=16, repr.plot.height=10)
p2 <- ggplot(biotur1, aes(x = as.integer(as.character(Time)), y = abs(Bio_rate_1sqm))) + 
  geom_line() +
  geom_point(size = 5) +
  facet_grid(~quadrat) +
  scale_x_continuous("Time elapsed (days)", breaks = c(1, 13, 14), limits = c(0, 16)) +
  scale_y_continuous("Absolute change of volume in \nthe sediment (cm3 per m2d-1)", breaks = seq(500, 3500, 500)) +
  theme(
    axis.text = element_text(size = 15),
    axis.title = element_text(size = 18),
    # legend.position = "none",
    panel.border = element_blank(),
    panel.background = element_blank(),
    panel.grid.minor = element_blank(),
    panel.grid.major = element_blank(),
    axis.line = element_line(color = "black", size = .5),
    # legend.key = element_rect(fill = "white", color = NA),
    axis.text.x = element_text(angle = 90, size = 14),
    strip.placement = "inside",
    strip.text = element_text(size = 18),
    strip.background = element_blank(),
    legend.text = element_text(size = 18),
    legend.title = element_text(size = 20))
p2

Aim 2: Study the potential process explaining bioturbation patterns

  • Question: Do the number of burrows explain the observed change in volume?
  • Question: Do the distribution of burrows in the quadrat is associated with areas of major change in volume?
  • Question: Do the burrows distribution correlates with the local geometry of the sediment surface? This is, burrows position explain the observed changes in verticality, surface variation and other geometric features.

#collapse-hide
## Load clean data
burrows <- readRDS(here::here('data/clean', 'burrows_master.rds'))

library(dplyr)

options(repr.plot.width=16, repr.plot.height=8)

burrows_tally <- burrows %>%
  group_by(ID) %>%
  tally() %>%
  mutate(quadrat = as.factor(as.character(str_remove_all(ID, "burrowmap_C01_")))) %>%
  mutate(quadrat = as.factor(as.character(stringr::str_sub(quadrat, 1, 1)))) %>%
  mutate(quadrat_time = as.factor(as.character(stringr::str_sub(ID, -3, -1)))) %>%
  mutate(quadrat_time = gsub("_", "", quadrat_time))

# burrows_tally$time <- rep(seq(1,3,1), 5)

biotur1$nburrows1 <- burrows_tally$n[match(interaction(burrows_tally$quadrat, burrows_tally$quadrat_time), 
                                          interaction(biotur1$quadrat, biotur1$time1))]
biotur1$nburrows2 <- burrows_tally$n[match(interaction(burrows_tally$quadrat, burrows_tally$quadrat_time), 
                                          interaction(biotur1$quadrat, biotur1$time2))]
biotur1$nburrows1 <- ifelse(is.na(biotur1$nburrows1), biotur1$nburrows2, biotur1$nburrows1)


biotur1$comparison <- paste0(biotur1$time1, "_vs_", biotur1$time2)



p3 <- ggplot() +
  geom_point(data = subset(biotur1, variable == "Removed.volume" & time2 == "D1"), 
            aes(x = nburrows1, y = log(abs(Bio_rate_1sqm)), colour = "col1"), size = 5) +
  geom_point(data = subset(biotur1, variable == "Removed.volume" & time1 == "D1" & time2 == "D14"), 
             aes(x = nburrows1, y = log(abs(Bio_rate_1sqm)), colour = "col2"), size = 5) +
  geom_point(data = subset(biotur1, variable == "Removed.volume" & time1 == "D0" & time2 == "D14"), 
             aes(x = nburrows1, y = log(abs(Bio_rate_1sqm)), colour = "col3"), size = 5) +
  scale_colour_manual("Time elapsed", labels = c("1 day", "13 day", "14 days"),
                      values = c("col1" = "#bf812d", "col2" = "#35978f", "col3" = "#000000")) + 
  labs(x = "Number of burrows", y = "Absolute change of volume in \nthe sediment ( log(cm3 per m2d-1) )") +
  theme(axis.text = element_text(size = 15),
    axis.title = element_text(size = 18),
    legend.position = c(0.8, 0.2),
    panel.border = element_blank(),
    panel.background = element_blank(),
    panel.grid.minor = element_blank(),
    panel.grid.major = element_blank(),
    axis.line = element_line(color = "black", size = .5),
    # legend.key = element_rect(fill = "white", color = NA),
    axis.text.x = element_text(size = 14),
    # strip.placement = "inside",
    # strip.text = element_text(size = 18),
    # strip.background = element_blank(),
    legend.text = element_text(size = 18),
    legend.title = element_text(size = 20))
  

p4 <- ggplot() +
    geom_point(data = subset(biotur1, variable == "Removed.volume" & time2 == "D1"), 
               aes(x = nburrows1, y = abs(Bio_rate_1sqm), colour = "col1"), size = 5) +
    geom_point(data = subset(biotur1, variable == "Removed.volume" & time1 == "D1" & time2 == "D14"),
               aes(x = nburrows1, y = abs(Bio_rate_1sqm), colour = "col2"), size = 5) +
    geom_point(data = subset(biotur1, variable == "Removed.volume" & time1 == "D0" & time2 == "D14"), 
               aes(x = nburrows1, y = abs(Bio_rate_1sqm), colour = "col3"), size = 5) +
    scale_colour_manual("Time elapsed", labels = c("1 day lapsed", "13 day lapsed", "14 days elapse"),
                        values = c("col1" = "#bf812d", "col2" = "#35978f", "col3" = "#000000")) + 
    labs(x = "Number of burrows", y = "Absolute change of volume in \nthe sediment (cm3 per m2d-1)") +
  theme(axis.text = element_text(size = 15),
    axis.title = element_text(size = 18),
    legend.position = "none",
    panel.border = element_blank(),
    panel.background = element_blank(),
    panel.grid.minor = element_blank(),
    panel.grid.major = element_blank(),
    axis.line = element_line(color = "black", size = .5),
    # legend.key = element_rect(fill = "white", color = NA),
    axis.text.x = element_text(size = 14),
    # strip.placement = "inside",
    # strip.text = element_text(size = 18),
    # strip.background = element_blank(),
    legend.text = element_text(size = 18),
    legend.title = element_text(size = 20))

plot_grid(p3, p4)

Discussion

Final remarks